{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The heat equation in 1D\n",
    "\n",
    "$$\n",
    "\\renewcommand{\\DdQq}[2]{{\\mathrm D}_{#1}{\\mathrm Q}_{#2}}\n",
    "\\renewcommand{\\drondt}{\\partial_t}\n",
    "\\renewcommand{\\drondx}{\\partial_x}\n",
    "\\renewcommand{\\drondtt}{\\partial_{tt}}\n",
    "\\renewcommand{\\drondxx}{\\partial_{xx}}\n",
    "\\renewcommand{\\drondyy}{\\partial_{yy}}\n",
    "\\renewcommand{\\dx}{\\Delta x}\n",
    "\\renewcommand{\\dt}{\\Delta t}\n",
    "\\renewcommand{\\grandO}{{\\mathcal O}}\n",
    "\\renewcommand{\\density}[2]{\\,f_{#1}^{#2}}\n",
    "\\renewcommand{\\fk}[1]{\\density{#1}{\\vphantom{\\star}}}\n",
    "\\renewcommand{\\fks}[1]{\\density{#1}{\\star}}\n",
    "\\renewcommand{\\moment}[2]{\\,m_{#1}^{#2}}\n",
    "\\renewcommand{\\mk}[1]{\\moment{#1}{\\vphantom{\\star}}}\n",
    "\\renewcommand{\\mke}[1]{\\moment{#1}{e}}\n",
    "\\renewcommand{\\mks}[1]{\\moment{#1}{\\star}}\n",
    "$$\n",
    "\n",
    "In this tutorial, we test a very classical lattice Boltzmann scheme $\\DdQq{1}{3}$ on the heat equation.\n",
    "\n",
    "The problem reads\n",
    "$$\n",
    "\\begin{gathered} \\drondt u = \\mu \\drondxx u, \\quad t>0, \\quad x\\in(0,1),\\\\ u(0) = u(1) = 0, \\end{gathered}\n",
    "$$\n",
    "where $\\mu$ is a constant scalar."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The scheme $\\DdQq{1}{3}$\n",
    "\n",
    "The numerical simulation of this equation by a lattice Boltzmann scheme consists in the approximatation of the solution on discret points of $(0,1)$ at discret instants.\n",
    "\n",
    "To simulate this system of equations, we use the $\\DdQq{1}{3}$ scheme given by\n",
    "\n",
    "* three velocities $v_0=0$, $v_1=1$, and $v_2=-1$, with associated distribution functions $\\fk{0}$, $\\fk{1}$, and $\\fk{2}$,\n",
    "* a space step $\\dx$ and a time step $\\dt$, the ration $\\lambda=\\dx/\\dt$ is called the scheme velocity,\n",
    "* three moments\n",
    "  $$ \\mk{0}=\\sum_{i=0}^{2} \\fk{i}, \\quad \\mk{1}= \\sum_{i=0}^{2} v_i \\fk{i}, \\quad \\mk{2}= \\frac{1}{2} \\sum_{i=0}^{2} v_i^2 \\fk{i},$$\n",
    "  \n",
    "  and their equilibrium values $\\mke{0}$, $\\mke{1}$, and $\\mke{2}$.\n",
    "* two relaxation parameters $s_1$ and $s_2$ lying in $[0,2]$.\n",
    "\n",
    "In order to use the formalism of the package pylbm, we introduce the three polynomials that define the moments: $P_0 = 1$, $P_1=X$, and $P_2=X^2/2$, such that\n",
    "$$ \\mk{k} = \\sum_{i=0}^2 P_k(v_i) \\fk{i}.$$\n",
    "\n",
    "The transformation $(\\fk{0}, \\fk{1}, \\fk{2})\\mapsto(\\mk{0},\\mk{1}, \\mk{2})$ is invertible if, and only if, the polynomials $(P_0,P_1,P_2)$ is a free set over the stencil of velocities.\n",
    "\n",
    "The lattice Boltzmann method consists to compute the distribution functions $\\fk{0}$, $\\fk{1}$, and $\\fk{2}$ in each point of the lattice $x$ and at each time $t^n=n\\dt$.\n",
    "A step of the scheme can be read as a splitting between the relaxation phase and the transport phase:\n",
    "\n",
    "* relaxation: \n",
    "$$\n",
    "    \\begin{aligned}\\mks{1}(t,x)&=(1-s_1)\\mk{1}(t,x)+s_1\\mke{1}(t,x),\\\\ \\mks{2}(t,x)&=(1-s_2)\\mk{2}(t,x)+s_2\\mke{2}(t,x).\\end{aligned}\n",
    "$$    \n",
    "\n",
    "* m2f: \n",
    "$$\n",
    "    \\begin{aligned}\\fks{0}(t,x)&\\;=\\mk{0}(t,x)-2\\mks{2}(t,x), \\\\ \\fks{1}(t,x)&\\;=\\mks{1}(t,x)/2+\\mks{2}(t,x), \\\\ \\fks{2}(t,x)&\\;=-\\mks{1}(t,x)/2+\\mks{2}(t,x).\\end{aligned}\n",
    "    $$\n",
    "\n",
    "* transport: \n",
    "$$\n",
    "    \\begin{aligned} \\fk{0}(t+\\dt, x)&\\;=\\fks{0}(t,x), \\\\ \\fk{1}(t+\\dt, x)&\\;=\\fks{1}(t,x-\\dx), \\\\ \\fk{2}(t+\\dt, x)&\\;=\\fks{2}(t,x+\\dx). \\end{aligned}\n",
    "    $$\n",
    "\n",
    "* f2m:\n",
    "$$\n",
    "    \\begin{aligned}\\mk{0}(t+\\dt,x)&\\;=\\fk{0}(t+\\dt,x)+\\fk{1}(t+\\dt,x)+\\fk{2}(t+\\dt,x), \\\\ \\mk{1}(t+\\dt,x)&\\;=\\fk{1}(t+\\dt,x)-\\fk{2}(t+\\dt,x), \\\\ \\mk{2}(t+\\dt,x)&\\;=\\tfrac{1}{2}\\fk{1}(t+\\dt,x)+\\tfrac{1}{2}\\fk{2}(t+\\dt,x).\\end{aligned}\n",
    "$$\n",
    "\n",
    "The moment of order $0$, $\\mk{0}$, being conserved during the relaxation phase, \n",
    "a diffusive scaling $\\dt=\\dx^2$, yields to the following equivalent equation\n",
    "$$\\drondt\\mk{0} = 2\\bigl(\\tfrac{1}{s_1}-\\tfrac{1}{2}\\bigr) \\drondxx\\mke{2} + \\grandO(\\dx^2),$$\n",
    "\n",
    "if $\\mke{1}=0$.\n",
    "In order to be consistent with the heat equation, the following choice is done:\n",
    "$$\\mke{2}=\\tfrac{1}{2}u, \\qquad s_1 = \\frac{2}{1+2\\mu}, \\qquad s_2=1.$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using pylbm\n",
    "\n",
    "pylbm uses Python dictionary to describe the simulation. In the following, we will build this dictionary step by step.\n",
    "\n",
    "### The geometry\n",
    "\n",
    "In pylbm, the geometry is defined by a box and a label for the boundaries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+----------------------+\n",
      "| Geometry information |\n",
      "+----------------------+\n",
      "    - spatial dimension: 1\n",
      "    - bounds of the box: [0. 1.]\n",
      "    - labels: [0, 0]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAADnCAYAAAC9roUQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAABdBJREFUeJzt3D9rVVkUxuF1RDQEK7k2ikyagAgOgtFWa8EmIAFTW9n45wsEv4DhaqWtwa9hbRKQQbCwiYU2E6YzRMG7pxgcjJOkSt5zJ/d5IEXOTrGKxY+dkHO71loBkHGs7wEAJonoAgSJLkCQ6AIEiS5A0PH9DgeDQZuZmQmNAnA0rK+vb7bWzux2tm90Z2Zmam1t7XCmAjiiuq77uNeZPy8ABIkuQJDoAgSJLkCQ6AIEiS5AkOgCBIkuQJDoAgSJLkCQ6AIEiS5AkOgCBIkuQJDoAgSJLkCQ6AIEiS5AkOgCBIkuQJDoAgSJLkCQ6AIEiS5AkOgCBIkuQJDoAgSJLkCQ6AIEiS5AkOgCBIkuQJDoAgSJLkCQ6AIEiS5AkOgCBIkuQJDoAgSJLkCQ6AIEiS5AkOgCBIkuQJDoAgSJLkCQ6AIEiS5AkOgCBIkuQJDoAgSJLkCQ6AIEiS5AkOgCBIkuQJDoAgSJLkCQ6AIEiS5AkOgCBIkuQJDoAgSJLkCQ6AIEiS5AkOgCBIkuQJDoAgSJLkCQ6AIEiS5AkOgCBIkuQJDoAgSJLkCQ6AIEiS5AkOgCBIkuQJDoAgSJLkCQ6AIEiS5AkOgCBIkuQJDoAgSJLkCQ6AIEiS5AkOgCBIkuQJDoAgSJLkCQ6AIEHe97gEmxsvJHffjw145nS0s3+hmGibe09HrH97Ozp2tx8fd+hpkwbrohvwYXxon9zHHTDXO7ZRz8vIe/3no5XG66AEGiCxAkugBBogsQJLoAQaILECS6AEGiCxDk5YiQ1683qso/ojN+fuwmGW66AEFuuiE3bsxUldeAGT9++8py0wUIEl2AINEFCBJdgCDRBQgSXYAg0QUIEl2AIC9HhP38j+helKAvXojoj5tuyOzs6b5HgD3Zz5yutbbn4dzcXFtbWwuOA/D/13XdemttbrczN12AINEFCBLdoNFoVE+ePKkLFy7U1NRUnT9/vh49elRfvnzpezQmgP0bD6Ib9ODBg3r48GFdvHixnj59Wrdv367hcFi3bt2q0WjU93gccfZvTLTW9vy6cuVK42C8e/eudV3X5ufndzwfDoetqtrKykpPkzEJ7F9WVa21Pbrqphvy6tWraq3V/fv3dzy/e/duTU9P18uXL3uajElg/8aH6Iasrq7WsWPH6tq1azueT01N1eXLl2t1dbWnyZgE9m98iG7I58+fazAY1MmTJ/9zdu7cudrc3Kxv3771MBmTwP6ND9EN2dra2nXhq/65bfz4GTgM9m98iG7I9PR0ff36ddez7e3tf38GDoP9Gx+iG3L27Nna3NzcdfE/ffpUg8GgTpw40cNkTAL7Nz5EN+Tq1as1Go3qzZs3O55vb2/X27dva25u19e04UDYv/EhuiELCwvVdV0tLy/veP7ixYva2tqqxcXFniZjEti/8eHzdEMuXbpU9+7dq2fPntX8/HzdvHmz3r9/X8PhsK5fv1537tzpe0SOMPs3Pny0Y9D3799reXm5nj9/XhsbGzUYDGphYaEeP35cp06d6ns8jjj7l7PfRzuKLsAB83m6AGNCdAGCRBcgSHQBgkQXIEh0AYJEFyBIdAGCRBcgSHQBgkQXIEh0AYJEFyBIdAGCRBcgSHQBgkQXIEh0AYJEFyBIdAGCRBcgSHQBgkQXIEh0AYJEFyBIdAGCRBcgSHQBgkQXIEh0AYJEFyBIdAGCRBcgSHQBgkQXIEh0AYJEFyBIdAGCRBcgSHQBgkQXIEh0AYJEFyBIdAGCRBcgSHQBgkQXIEh0AYJEFyBIdAGCRBcgSHQBgkQXIEh0AYJEFyBIdAGCRBcgSHQBgkQXIEh0AYJEFyBIdAGCRBcgSHQBgkQXIEh0AYJEFyBIdAGCRBcgSHQBgkQXIEh0AYJEFyBIdAGCRBcgSHQBgkQXIEh0AYJEFyBIdAGCutba3odd92dVfcyNA3Ak/NZaO7Pbwb7RBeBg+fMCQJDoAgSJLkCQ6AIEiS5A0N9fpYJgC2jnKgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import pylbm\n",
    "import numpy as np\n",
    "\n",
    "xmin, xmax = 0., 1.\n",
    "dico_geom = {\n",
    "    'box': {'x': [xmin, xmax], 'label':0},\n",
    "}\n",
    "geom = pylbm.Geometry(dico_geom)\n",
    "print(geom)\n",
    "geom.visualize(viewlabel=True);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The stencil\n",
    "\n",
    "pylbm provides a class stencil that is used to define the discret velocities of the scheme. In this example, the stencil is composed by the velocities $v_0=0$, $v_1=1$ and $v_2=-1$ numbered by $[0,1,2]$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+---------------------+\n",
      "| Stencil information |\n",
      "+---------------------+\n",
      "    - spatial dimension: 1\n",
      "    - minimal velocity in each direction: [-1]\n",
      "    - maximal velocity in each direction: [1]\n",
      "    - information for each elementary stencil:\n",
      "        stencil 0\n",
      "            - number of velocities: 3\n",
      "            - velocities\n",
      "                (0: 0)\n",
      "                (1: 1)\n",
      "                (2: -1)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATAAAAC6CAYAAAAppxH5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAACs1JREFUeJzt3H+slfV9wPH3h18We8HrDwqUH0E3KkUn0jGm6B+mXat1y8QuS2qzpUzn0qQ2pXFZu3Q/XJc16VarZpiyalmdmLYma9Vs1f7IahqsWvk9b0HiFkFEoBeHghaF62d/nIO52+64F85tnvu59/1KCOc55znP/eSbwzvPec65RGYiSRWNa3oASTpVBkxSWQZMUlkGTFJZBkxSWQZMUlkGTCNCRMyNiMMRMb69/WhE/GHTc2lkM2AiIi6PiB9HxMsR8VJEPBYRv9Z+bEVErPtFz5CZuzKzKzP7hrJ/RHwqIva2Z14TEaf9omfUyGPAxriImAr8C/D3wFnALOCvgNebnOtEIuJK4DPA+4B5wHm0ZtYYY8D0LoDM/Hpm9mXmzzPze5m5NSLeDawGLm2/vTsIEBGnRcQXI2JXROyLiNURMbn92BURsTsibo6I/RHxYkT8wfEfFhGTI+LWiNjZPnta175vXkRkREwYwswfBb6amT2Z+V/AXwMrhnldVIAB0w6gLyLuiYgPRsSZxx/IzG3Ax4DH22/vutsPfYFW+C4GfpnWWdtf9DvmDOCM9v03AHf2O+4XgV8FltE64/sT4M2TnPkCYEu/7S3A9Ig4+ySPo+IM2BiXma8AlwMJ3AX8LCIeiojpA+0fEQHcCHwqM1/KzEPA54EP99vtKPC5zDyamd8BDgPnR8Q44Hrgk5n5QvuM78eZebJvV7uAl/ttH7895SSPo+KGcrquUa59prUCICIWAGuB24HrBth9GnA6sKHVMgACGN9vnwOZeazf9mu0onMO8DbgPzoc+TAwtd/28duHOjyuivEMTP9DZm4HvgZcePyu/7VLL/Bz4ILM7G7/OSMzu4Zw+F7gCPBLHY7ZAyzqt70I2JeZBzo8rooxYGNcRCxoX3Cf3d6eQ+vM64n2LvuA2RExCSAz36T1VvO2iHhH+zmz2p8MnlD7uWuAL0XEOyNifERcegpfgfgn4IaIWNi+tvZntKKrMcaA6RDw68CTEfEqrXA9DdzcfvzfaJ3x7I2I3vZ9nwaeBZ6IiFeAHwDnD/Hn/THw78BTwEu0PhA4qddhZj4C/C3wQ2Bn+89fnswxNDqE/6GhpKo8A5NUlgGTVJYBk1SWAZNUVkdfZD377LPz3HPPHa5ZRoW+vj7Gjx8/+I5jjOsyMNdlYBs2bOjNzGmD7ddRwObOncv69es7OcSoc/DgQbq7uwffcYxxXQbmugwsInYOZT/fQkoqy4BJKsuASSrLgEkqy4BJKsuASSrLgEkqy4BJKsuASSrLgEkqy4BJKsuASSrLgEkqy4BJKsuASSrLgEkqy4BJKsuASSrLgEkqy4BJKsuASSrLgEkqy4BJKsuASSrLgEkqy4BJKsuASSrLgEkqy4BJKsuASSrLgEkqy4BJKsuASSrLgEkqy4BJKsuASSrLgEkqy4BJKsuASSrLgEkqy4BJKsuASSrLgEkqy4BJKsuASSrLgEkqy4BJKsuASSrLgEkqy4BJKsuASSrLgEkqy4BJKsuASSrLgEkqy4BJKsuASSrLgEkqy4BJKsuASSrLgEkqy4BJKsuASSrLgA2DG298iIUL76Sr6/Ocd94/cPXV99HTs7/psUaEI0eO8YlPfIf587/C5Ml/w2WXreHJJ3c3PVajbr/9CS666MuMH/85zjzzDm655dGmRyrLgA2Du+/exNSpp3HddRcyZcokHn74Wa68ci1HjhxrerTGrVz5CKtWPcW0aaezfPkCHn/8ed7//nvp7X2t6dEas2HDi5x11mTmzJna9CjlTWh6gNHgsceuZ9myOQBs3bqLRYv+kRdeOMRPf/oz3vOemQ1P15z9+19lzZpNjBsXPPjgh5g//51MmDCOtWu3smrVT7jlliuaHrER9957LQDLl3+DnTtfbnia2jwDGwbH4wXwxht9AIwbF8yc2dXUSCNCT89+jh59k7lzz2DatNMBWLKkFfTNm/c2OZpGiY7OwOLYYdj36DCNUt/hV/v4+B9tAuDmj81i5rgNsK/hoRq0b0frOmDX244y4cA6eL2Lt/e9CMDe5/f42nm9t/X34edci1PkGdgw6T1wlPf+zhZ+sulVbvy9mXzhz89reqTGTT9nItAK+3HHb894x6RGZtLo0tEZWE7ogulXDNMode3ceZAPXLuWHTsOsXLlEm677TebHmlEWHj5YSZOfJpdL7zOnr7FvGv6LJ7a/i0AFi29wNfOaXuBA9A1z7U4RV7EHwbLlq1hz55DzJ17BkeOHGPlykcA+MhHfoWlS2c1PF1zpk/vYsWKi7nrro1cc823WbRoBvff30NX1yRuumlp0+M15u67N7Ju3S42bmy9nX7gge0899xBli9fwPLlCxqerhYDNgz27DkEwK5dL7N69ea37r/44hljOmAAd9xxFRMnjuOb3+zhgQe2c8kls7n11g8wbdrbmx6tMevW7eKee7a8tb1lyz62bNnHvHndBuwkRWae8pMXL16cmzZtGsZx6jt48CDd3d1NjzHiuC4Dc10GFhEbMnPJYPt5EV9SWQZMUlkGTFJZBkxSWQZMUlkGTFJZBkxSWQZMUlkGTFJZBkxSWQZMUlkGTFJZBkxSWQZMUlkGTFJZBkxSWQZMUlkGTFJZBkxSWQZMUlkGTFJZBkxSWQZMUlkGTFJZBkxSWQZMUlkGTFJZBkxSWQZMUlkGTFJZBkxSWQZMUlkGTFJZBkxSWQZMUlkGTFJZBkxSWQZMUlkGTFJZBkxSWQZMUlkGTFJZBkxSWQZMUlkGTFJZBkxSWQZMUlkGTFJZBkxSWQZMUlkGTFJZBkxSWQZMUlkGTFJZBkxSWQZMUlkGTFJZBkxSWQZMUlkGTFJZBkxSWQZMUlkGTFJZBkxSWQZMUlmRmaf+5IhDwDPDN86ocA7Q2/QQI5DrMjDXZWDnZ+aUwXaa0OEPeSYzl3R4jFElIta7Jv+X6zIw12VgEbF+KPv5FlJSWQZMUlmdBuwrwzLF6OKaDMx1GZjrMrAhrUtHF/ElqUm+hZRUlgGTVFZHAYuIv4uI7RGxNSK+HRHdwzVYZRHxuxHRExFvRsSY/4g8Iq6KiGci4tmI+EzT84wEEbEmIvZHxNNNzzJSRMSciPhhRGxr//v55GDP6fQM7PvAhZl5EbAD+NMOjzdaPA18CPhR04M0LSLGA3cCHwQWAtdFxMJmpxoRvgZc1fQQI8wx4ObMfDdwCfDxwV4rHQUsM7+Xmcfam08Aszs53miRmdsy099QaFkKPJuZ/5mZbwDfAK5peKbGZeaPgJeanmMkycwXM3Nj+/YhYBsw60TPGc5rYNcDDw/j8TQ6zAKe77e9m0FelFJEzAMWA0+eaL9Bf5UoIn4AzBjgoc9m5oPtfT5L6/TvvpMdtKqhrIsAiAHu87s7+n9FRBfwz8DKzHzlRPsOGrDM/I1BfthHgd8C3pdj6Etlg62L3rIbmNNvezawp6FZNMJFxERa8bovM7812P6dfgp5FfBp4Lcz87VOjqVR6ylgfkScGxGTgA8DDzU8k0agiAjgq8C2zPzSUJ7T6TWwVcAU4PsRsTkiVnd4vFEhIq6NiN3ApcC/RsR3m56pKe0PeW4Cvkvrouz9mdnT7FTNi4ivA48D50fE7oi4oemZRoDLgN8H3tvuyeaIuPpET/BXiSSV5TfxJZVlwCSVZcAklWXAJJVlwCSVZcAklWXAJJX131h5o7gzTnV4AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "dico_sten = {\n",
    "    'dim': 1,\n",
    "    'schemes':[{'velocities': list(range(3))}],\n",
    "}\n",
    "sten = pylbm.Stencil(dico_sten)\n",
    "print(sten)\n",
    "sten.visualize();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The domain\n",
    "\n",
    "In order to build the domain of the simulation, the dictionary should contain the space step $\\dx$ and the stencils of the velocities (one for each scheme). \n",
    "\n",
    "We construct a domain with $N=10$ points in space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+--------------------+\n",
      "| Domain information |\n",
      "+--------------------+\n",
      "    - spatial dimension: 1\n",
      "    - space step: 0.1\n",
      "    - with halo:\n",
      "        bounds of the box: [-0.05] x [1.05]\n",
      "        number of points: [12]\n",
      "    - without halo:\n",
      "        bounds of the box: [0.05] x [0.95]\n",
      "        number of points: [10]\n",
      "        \n",
      "    +----------------------+\n",
      "    | Geometry information |\n",
      "    +----------------------+\n",
      "        - spatial dimension: 1\n",
      "        - bounds of the box: [0. 1.]\n",
      "        - labels: [0, 0]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWQAAACeCAYAAAAFZXDrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAACrBJREFUeJzt3G2MpWV9x/Hvb91QQllAZuShCkxpIakhLZoVrWlLjdogWvWFqbRF14Y0ARPfmCal0RdG08Q2QZNWE9tSrQ9BqU1r1z4kooXagKsOQRclkSdZ3QILu+oWXVax8++L+55wOszOObN7Hq6z+/0kk7nPOddc53fu3Oc397nOmUlVIUmavS2zDiBJ6ljIktQIC1mSGmEhS1IjLGRJaoSFLEmNsJB1QkpyfpIfJnnWrLNIqyxkTVWSh5I8meSJJD9IckeSa5NM9Visqu9U1alV9b/TvF9pIxayZuG3q2obcAHwXuCPgb+dbSRp9ixkzUxVHayqncAbgR1JLklyepKPJXk8yZ4k71w9e07yliS3J3l/f3b9YJKX9td/N8ljSXaszp/k1UnuSvI//e3vGrhtKUkl2dpfvi3Je/r5n0jyuSSLU94lOsFZyJq5qvoKsBf4deAvgdOBC4HLgTcDfzAw/MXAbmABuAn4FPAi4BeBq4EPJDm1H/uj/ufPAF4NXJfk9RtE+b3+vs4CTgL+aAwPTxqZhaxWPAycSXe2/CdV9URVPQTcALxpYNy3q+oj/drvzcB5wLur6sdV9TngJ3TlTFXdVlV3V9VKVe0GPklX8kfykaq6t6qeBP4euHTMj1HakIWsVjwX2Ep3Zrpn4Po9/W2r9g1sPwlQVWuvOxUgyYuT3NovfxwErgU2WoZ4dGD70Oo80rRYyJq5JC+iK93PAE/Rvdm36nzgv49y6puAncB5VXU68CEgxxBVmigLWTOT5LQkr6FbB/5EVX2dbqngT5NsS3IB8HbgE0d5F9uA71XV4SSX0a0RS83aOusAOiF9NslPgRXgHuB9dGevAG+je2PvQeAw8DfAh4/yft4K3JDkA8B/0pX9GceQW5qo+A/qJakNLllIUiMsZElqhIUsSY2wkCWpEZv6lMXi4mItLS1NKIokHZ/uvPPO/VX1nGHjNlXIS0tLLC8vH30qSToBJdkzfJRLFpLUDAtZkhphIUtSIyxkSWqEhSxJjbCQJakRFrIkNcJClqRGWMiS1AgLWZIaYSFLUiMsZElqhIUsSY2wkCWpERayJDXCQpakRljIktQIC1mSGmEhS1IjLGRJaoSFLEmNsJAlqREWsiQ1wkKWpEZYyJLUCAtZkhphIUtSIyxkSWqEhSxJjbCQJakRFrIkNcJClqRGWMiS1AgLWZIaYSFLUiMsZElqhIUsSY2wkCWpERayJDXCQpakRljIktQIC1mSGmEhS1IjLGRJaoSFLEmNsJAlqREWsiQ1wkKWpEZYyJLUCAtZkhphIUtSIyxkSWrE1lkHmIi7rofDjz59+eRz4AXvHdv0+/cf4sCBQywsnMLi4ikn7LyTnHve5p3k3PM278TmnvDzugXHZyEffhR+dunpyz96aGxT799/iFtueYAtW8LKSvHKV/7CWA64eZt3knPP27yTnHve5p3o3BN8XrfCJYtNOnDgEFu2hLPPPpUtW8KBA4dOyHknOfe8zTvJuedt3knPfbyzkDdpYeEUVlaKfft+yMpKsbAwnrOKeZt3knPP27yTnHve5p303Me7VNXIg7dv317Ly8sTjDMmriFPZd5Jzj1v805y7nmbd2Jzz/EacpI7q2r70HHHZSFLUkNGLWSXLCSpERayJDXCQpakRljIktQIC1mSGmEhS1IjLGRJaoSFLEmNsJAlqREWsiQ1wkKWpEZYyJLUCAtZkhphIUtSIyxkSWqEhSxJjbCQJakRFrIkNcJClqRGWMiS1AgLWZIaYSFLUiMsZElqhIUsSY2wkCWpERayJDXCQpakRljIktQIC1mSGmEhS1IjLGRJaoSFLEmNsJAlqREWsiQ1wkKWpEZYyJLUCAtZkhphIUtSIyxkSWqEhSxJjbCQJakRFrIkNcJClqRGWMiS1AgLWZIaYSFLUiMsZElqhIUsSY2wkCWpERayJDXCQpakRljIktSIVNXog5PHgT0jDl8E9h9NqBmat8zzlhfmL/O85QUzT8Nm815QVc8ZNmhThbwZSZaravtEJp+Qecs8b3lh/jLPW14w8zRMKq9LFpLUCAtZkhoxyUL+6wnOPSnzlnne8sL8ZZ63vGDmaZhI3omtIUuSNsclC0lqhIUsSY0YWyEnOTPJLUnu678/e50xlyb5UpJvJtmd5I3juv9N5LwiybeS3J/k+nVu/5kkN/e3fznJ0rQzrpNpWOa3J7mn36dfSHLBLHKuybRh5oFxb0hSSWb6kadR8ib5nX4/fzPJTdPOuE6eYcfF+UluTXJXf2xcOYucA3k+nOSxJN84wu1J8hf949md5IXTzrgmz7C8v9/n3J3kjiS/csx3WlVj+QL+HLi+374e+LN1xlwMXNRv/xzwCHDGuDKMkPFZwAPAhcBJwNeB568Z81bgQ/32VcDN08p3DJlfBpzSb183D5n7cduALwK7gO0t5wUuAu4Cnt1fPqv1fUz3xtN1/fbzgYdmnPk3gBcC3zjC7VcC/w4EeAnw5cbzvnTgeHjVOPKOc8nidcBH++2PAq9fO6Cq7q2q+/rth4HHgKF/vTJGlwH3V9WDVfUT4FN0uQcNPo5/AF6eJFPMuNbQzFV1a1Ud6i/uAp435YxrjbKfAd5D94v88DTDrWOUvH8IfLCqvg9QVY9NOeNao2Qu4LR++3Tg4Snme4aq+iLwvQ2GvA74WHV2AWckOXc66Z5pWN6qumP1eGBMz7txFvLZVfUIQP/9rI0GJ7mM7jf7A2PMMMxzge8OXN7bX7fumKr6KXAQWJhKuvWNknnQNXRnGbM0NHOSFwDnVdW/TDPYEYyyjy8GLk5ye5JdSa6YWrr1jZL5XcDVSfYC/wa8bTrRjtpmj/WWjOV5t3Uzg5N8HjhnnZvescl5zgU+DuyoqpXN/OwxWu9Md+3n/kYZM00j50lyNbAduHyiiYbbMHOSLcD7gbdMK9AQo+zjrXTLFr9Jdyb0X0kuqaofTDjbkYyS+XeBv6uqG5L8KvDxPvM0n3Ob0dpzbyRJXkZXyL92rHNtqpCr6hUbhNqX5NyqeqQv3HVf0iU5DfhX4J39y5Jp2gucN3D5eTzzZdzqmL1JttK91NvoZdakjZKZJK+g+8V4eVX9eErZjmRY5m3AJcBt/WrQOcDOJK+tquWppXzaqMfFrqp6Cvh2km/RFfRXpxPxGUbJfA1wBUBVfSnJyXT/FGfWyy1HMtKx3pIkvwzcCLyqqg4c63zjXLLYCezot3cA/7x2QJKTgH+iWyf69Bjve1RfBS5K8vN9lqvocg8afBxvAP6j+lX7GRmauX/5/1fAaxtY24QhmavqYFUtVtVSVS3Rrb/NqoxhtOPiM3RvnpJkkW4J48Gppvz/Rsn8HeDlAEl+CTgZeHyqKTdnJ/Dm/tMWLwEOri6DtijJ+cA/Am+qqnvHMukY35FcAL4A3Nd/P7O/fjtwY799NfAU8LWBr0un/M7plcC9dGvX7+ivezddIUB30H4auB/4CnDhNPMdZebPA/sG9unO1jOvGXsbM/yUxYj7OMD7gHuAu4GrWt/HdJ+suJ3uExhfA35rxnk/SffJqqfozoavAa4Frh3Yxx/sH8/dDRwTw/LeCHx/4Hm3fKz36Z9OS1Ij/Es9SWqEhSxJjbCQJakRFrIkNcJClqRGWMiS1AgLWZIa8X8IPZ2XK4bmnAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "N = 10\n",
    "dx = (xmax-xmin)/N\n",
    "dico_dom = {\n",
    "    'box': {'x': [xmin, xmax], 'label':0},\n",
    "    'space_step': dx,\n",
    "    'schemes': [\n",
    "        {\n",
    "            'velocities': list(range(3)),\n",
    "        }\n",
    "    ],\n",
    "}\n",
    "dom = pylbm.Domain(dico_dom)\n",
    "print(dom)\n",
    "dom.visualize();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The scheme\n",
    "\n",
    "In pylbm, a simulation can be performed by using several coupled schemes. In this example, a single scheme is used and defined through a list of one single dictionary. This dictionary should contain:\n",
    "\n",
    "* 'velocities': a list of the velocities\n",
    "* 'conserved_moments': a list of the conserved moments as sympy variables\n",
    "* 'polynomials': a list of the polynomials that define the moments\n",
    "* 'equilibrium': a list of the equilibrium value of all the moments\n",
    "* 'relaxation_parameters': a list of the relaxation parameters ($0$ for the conserved moments)\n",
    "* 'init': a dictionary to initialize the conserved moments\n",
    "\n",
    "(see the documentation for more details)\n",
    "\n",
    "The scheme velocity could be taken to $1/\\dx$ and the inital value of $u$ to\n",
    "\n",
    "$$ u(t=0,x) = \\sin(\\pi x).$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+--------------------+\n",
      "| Scheme information |\n",
      "+--------------------+\n",
      "    - spatial dimension: 1\n",
      "    - number of schemes: 1\n",
      "    - number of velocities: 3\n",
      "    - conserved moments: [u]\n",
      "        \n",
      "    +----------+\n",
      "    | Scheme 0 |\n",
      "    +----------+\n",
      "        - velocities\n",
      "            (0: 0)\n",
      "            (1: 1)\n",
      "            (2: -1)\n",
      "\n",
      "        - polynomials\n",
      "                    \n",
      "            ⎡1 ⎤\n",
      "            ⎢  ⎥\n",
      "            ⎢X ⎥\n",
      "            ⎢  ⎥\n",
      "            ⎢ 2⎥\n",
      "            ⎢X ⎥\n",
      "            ⎢──⎥\n",
      "            ⎣2 ⎦\n",
      "\n",
      "        - equilibria\n",
      "                    \n",
      "            ⎡  u  ⎤\n",
      "            ⎢     ⎥\n",
      "            ⎢ 0.0 ⎥\n",
      "            ⎢     ⎥\n",
      "            ⎣0.5⋅u⎦\n",
      "\n",
      "        - relaxation parameters\n",
      "                    \n",
      "            ⎡       0.0       ⎤\n",
      "            ⎢                 ⎥\n",
      "            ⎢0.666666666666667⎥\n",
      "            ⎢                 ⎥\n",
      "            ⎣       1.0       ⎦\n",
      "\n",
      "    - moments matrices\n",
      "                \n",
      "        ⎡1  1    1 ⎤\n",
      "        ⎢          ⎥\n",
      "        ⎢0  10  -10⎥\n",
      "        ⎢          ⎥\n",
      "        ⎣0  50  50 ⎦\n",
      "\n",
      "    - inverse of moments matrices\n",
      "                \n",
      "        ⎡1    0    -1/50⎤\n",
      "        ⎢               ⎥\n",
      "        ⎢0  1/20   1/100⎥\n",
      "        ⎢               ⎥\n",
      "        ⎣0  -1/20  1/100⎦\n",
      "\n",
      "    \n"
     ]
    }
   ],
   "source": [
    "import sympy as sp\n",
    "\n",
    "def solution(x, t):\n",
    "    return np.sin(np.pi*x)*np.exp(-np.pi**2*mu*t)\n",
    "\n",
    "# parameters\n",
    "mu = 1.\n",
    "la = 1./dx\n",
    "s1 = 2./(1+2*mu)\n",
    "s2 = 1.\n",
    "u, X = sp.symbols('u, X')\n",
    "\n",
    "dico_sch = {\n",
    "    'dim': 1,\n",
    "    'scheme_velocity': la,\n",
    "    'schemes':[\n",
    "        {\n",
    "            'velocities': list(range(3)),\n",
    "            'conserved_moments': u,\n",
    "            'polynomials': [1, X, X**2/2],\n",
    "            'equilibrium': [u, 0., .5*u],\n",
    "            'relaxation_parameters': [0., s1, s2],\n",
    "        }\n",
    "    ],\n",
    "}\n",
    "\n",
    "sch = pylbm.Scheme(dico_sch)\n",
    "print(sch)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The simulation\n",
    "\n",
    "A simulation is built by defining a correct dictionary.\n",
    "\n",
    "We combine the previous dictionaries to build a simulation. In order to impose the homogeneous Dirichlet conditions in $x=0$ and $x=1$, the dictionary should contain the key 'boundary_conditions' (we use pylbm.bc.Anti_bounce_back function)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+------------------------+\n",
      "| Simulation information |\n",
      "+------------------------+\n",
      "        \n",
      "    +--------------------+\n",
      "    | Domain information |\n",
      "    +--------------------+\n",
      "        - spatial dimension: 1\n",
      "        - space step: 0.1\n",
      "        - with halo:\n",
      "            bounds of the box: [-0.05] x [1.05]\n",
      "            number of points: [12]\n",
      "        - without halo:\n",
      "            bounds of the box: [0.05] x [0.95]\n",
      "            number of points: [10]\n",
      "            \n",
      "        +----------------------+\n",
      "        | Geometry information |\n",
      "        +----------------------+\n",
      "            - spatial dimension: 1\n",
      "            - bounds of the box: [0. 1.]\n",
      "            - labels: [0, 0]\n",
      "        \n",
      "    +--------------------+\n",
      "    | Scheme information |\n",
      "    +--------------------+\n",
      "        - spatial dimension: 1\n",
      "        - number of schemes: 1\n",
      "        - number of velocities: 3\n",
      "        - conserved moments: [u]\n",
      "            \n",
      "        +----------+\n",
      "        | Scheme 0 |\n",
      "        +----------+\n",
      "            - velocities\n",
      "                (0: 0)\n",
      "                (1: 1)\n",
      "                (2: -1)\n",
      "\n",
      "            - polynomials\n",
      "                        \n",
      "                ⎡1 ⎤\n",
      "                ⎢  ⎥\n",
      "                ⎢X ⎥\n",
      "                ⎢  ⎥\n",
      "                ⎢ 2⎥\n",
      "                ⎢X ⎥\n",
      "                ⎢──⎥\n",
      "                ⎣2 ⎦\n",
      "\n",
      "            - equilibria\n",
      "                        \n",
      "                ⎡  u  ⎤\n",
      "                ⎢     ⎥\n",
      "                ⎢ 0.0 ⎥\n",
      "                ⎢     ⎥\n",
      "                ⎣0.5⋅u⎦\n",
      "\n",
      "            - relaxation parameters\n",
      "                        \n",
      "                ⎡       0.0       ⎤\n",
      "                ⎢                 ⎥\n",
      "                ⎢0.666666666666667⎥\n",
      "                ⎢                 ⎥\n",
      "                ⎣       1.0       ⎦\n",
      "\n",
      "        - moments matrices\n",
      "                    \n",
      "            ⎡1  1    1 ⎤\n",
      "            ⎢          ⎥\n",
      "            ⎢0  10  -10⎥\n",
      "            ⎢          ⎥\n",
      "            ⎣0  50  50 ⎦\n",
      "\n",
      "        - inverse of moments matrices\n",
      "                    \n",
      "            ⎡1    0    -1/50⎤\n",
      "            ⎢               ⎥\n",
      "            ⎢0  1/20   1/100⎥\n",
      "            ⎢               ⎥\n",
      "            ⎣0  -1/20  1/100⎦\n",
      "\n",
      "        \n"
     ]
    }
   ],
   "source": [
    "dico = {\n",
    "    'box': {'x':[xmin, xmax], 'label':0},\n",
    "    'space_step': dx,\n",
    "    'scheme_velocity': la,\n",
    "    'schemes':[\n",
    "        {\n",
    "            'velocities': list(range(3)),\n",
    "            'conserved_moments': u,\n",
    "            'polynomials': [1, X, X**2/2],\n",
    "            'equilibrium': [u, 0., .5*u],\n",
    "            'relaxation_parameters': [0., s1, s2],\n",
    "        }\n",
    "    ],\n",
    "    'init': {u:(solution,(0.,))},\n",
    "    'boundary_conditions': {\n",
    "        0: {'method': {0: pylbm.bc.AntiBounceBack,}},\n",
    "    },\n",
    "    'generator': 'numpy'\n",
    "}\n",
    "\n",
    "sol = pylbm.Simulation(dico)\n",
    "print(sol)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run a simulation\n",
    "\n",
    "Once the simulation is initialized, one time step can be performed by using the function one_time_step.\n",
    "\n",
    "We compute the solution of the heat equation at $t=0.1$. And, on the same graphic, we plot the initial condition, the exact solution and the numerical solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEICAYAAACktLTqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3XlUVdX7x/H3ZkZABAcUZ3OeQJyH1NJKy6nS8qvNWmqT81RZmjmkpmZmWmpWmkNqRuaUCqg4ggoqqOCMqKACioIM9/n9AfpDU0EEzgX2a627Fveefc/5cMvnbvbZ52wlImiapmkFi4XRATRN07Scp4u7pmlaAaSLu6ZpWgGki7umaVoBpIu7pmlaAaSLu6ZpWgGki7umZaCU+kQpNd/oHJr2uHRx11BKnVZKtb/ntbeUUjtyaP+ilKqaE/vKSUqptkqpiIyvichEEembC8fK9uep0nytlLqS/piilFIPaFtGKeWtlIpM/9wr3bPdVim1UCl1TSl1USk15J7t7ZRSR5VSN5VSPkqpitnJrBlPF3dNM3/vAd0AD6A+0Ano94C2JmAD8PIDto8FqgEVgaeAEUqpDgBKqRLAamAM4AoEAMtz5DfQ8p6I6EchfwCngfb3vPYWsCPDc3dgFRANnAI+zrCtCbALiAUuALMBm/Rt2wABbgDxwKsPyPAOEArEABuBihm2PQMcBeLS9+0H9E3fNhZYnKFtpfTjWaU/fzt9v9eBk0C/9NcdgATSimF8+sP9PvvrAhxJ/918gVr3fG7DgOD0bMsBu/v8brWARCA1/Tixj/jfZyfwXobnfYDdmbzHKv1zqHTP6+eBZzM8Hw8sS//5PWBnhm23P6OaRv8/qh+P/tA9dy1TSikL4G8gCCgLtAMGKaWeS2+SCgwGSgDN07e/DyAirdPbeIiIo4j8pyeolOoGfAK8BJQEtgNL07eVIO1L5bP0/Z8AWj5C/CjSerpFSSv0M5RSXiJyA+gIRKbnchSRyHtyVU/PMSg91zrgb6WUTYZmrwAdgMqk9arfujeAiIQC/YFd6ccplr7/UUqp2Ac9MuyiDmmf/W1B6a89EqWUC2lfYA/a113HSf+MTmTnWJrxdHHXbltzT2GZk2FbY6CkiHwpIkkichL4CegJICKBIrJbRFJE5DQwD2jzCMfuB0wSkVARSQEmAp7p473PAyEislJEkoGZwMWs7lhE/hGRE5LGD9gEPJnFt78K/CMi/6YfexpgD7TI0GaWiESKyFXSvgA9HyHbZBEp9qBHhqaOpP1lcFsc4PigcfeHcMzw/oz7cnrAce7druUjurhrt3W7p7C8n2FbRcD9nuL/CeAGaT1cpdTa9BN010grziUe4dgVgW8z7PsqoEj7K8EdOHe7oYhIxueZUUp1VErtVkpdTd/384+QzR04k+HYpvRjl83QJuMXzU3+v4DmpHjS/vK4rSgQn/5ZPOp+br8/476uP+A4927X8hFd3LWsOAecuqdn6SQiz6dv/4G0MfFqIlKUtML/KL3Kc6SNhWfcv72I7CRtDL/87YbpvdXyGd57AyiS4XnpDG1tSRvSmQa4pX9prcuQLbPiGEnaF8+9xz7/CL/bbf85Vvq0y/gHPTI0PULaydTbPNJfe7QAIjGkfZ4P2tddx1FKOQBPZOdYmvF0cdeyYi9wTSk1Uillr5SyVErVVUo1Tt/uBFwD4pVSNYEB97z/ElDlIfufC4xWStUBUEo5K6V6pG/7B6ijlHpJKWUFfEyGAg4cBForpSoopZyB0Rm22QC2pJ0ETlFKdQSevSdX8fT33c8K4IX06YHWwFDgFmknOB/VJaBcxvF6SZt26figR4b3/goMUUqVVUq5p+dY9KADKaXsSPu9AWzTn2fc12dKKZf0/1bvZtjXn0BdpdTL6e/5HAgWkaPZ+H01g+nirmVKRFKBzqSNJ58CLgPzgdtFcRjQi7Q/33/iv9PnxgK/pA+7vHKf/f8JfA0sSx/WOUzayU5E5DLQA5gMXCFtGp9/hvf+m368YCAQWJth23XSvgxWkDYLpxfgnWH7UdJOmJ5Mz+Z+T65jwGvAd+m/c2egs4gkZfKR3c9W0nrAF5VSlx/xvfNIG88/RNpn80/6awCk9/QznkdI4P+HYI6mP7/tC9JOkp4hbdbRVBHZACAi0aRNoZxA2ufVlPTzKlr+ox592E7TjKWU8iVtuqK+klTTHkD33DVN0wogXdw1TdMKID0so2maVgDpnrumaVoBZGXUgUuUKCGVKlUy6vCapmn5UmBg4GURKZlZO8OKe6VKlQgICDDq8JqmafmSUupM5q30sIymaVqBpIu7pmlaAaSLu6ZpWgFk2Ji7pmlaZpKTk4mIiCAxMdHoKHnOzs6OcuXKYW1tna336+KuaZrZioiIwMnJiUqVKvHot6/Pv0SEK1euEBERQeXKlbO1j0yHZdIX041SSh1+wHallJqllApXSgUrpbyylUTTNO0eiYmJFC9evFAVdgClFMWLF3+sv1iyMua+iLRlxB6kI2l36qtG2hqMP2Q7jaZp2j0KW2G/7XF/70yHZURkm1Kq0kOadAV+TV8VZrdSqphSqoyIXHisZJqWi27dusXp06eJjIwkMjKSmJgYrl+/TkLC/98d197eHicnJ1xdXXF3d6ds2bJUqlQp22OgmpaXcmLMvSx3L3sWkf7af4q7Uuo90nr3VKhQIQcOrWmZu3btGv7+/uzbt4/AwECOHDnCqVOnMJlM922vlOJB91yytLSkSpUq1K1bl4YNG9K4cWNatGiBo2NurK6nmYMWLVqwc+fD12fp27cvQ4YMoXbt2kycOJFPPvnkkd7v6OhIfHz8Q9s8qizdOCy9575WROreZ9s/pC1uvCP9+RZghIgEPmyfjRo1En2FqpYbTCYT+/btY+3atWzcuJHAwEBMJhNKKapXr46Hhwc1atSgatWqlC9fnjJlyuDq6krRokWxtbW9U9wTExO5fv06ly9f5sKFC0RERBAeHs7Ro0cJCgoiLCwMACsrKxo1akSHDh3o1KkTXl5ehXYoIaeFhoZSq1Yto2M8kuwU6ge9536/v1IqUEQaZbbPnOi5R3D3mpblSFt7UtPyVGhoKL/88gu//fYbkZGRWFpa0rx5cz799FPatGlDkyZNcHJyytK+lFLY29tjb29PqVKlqF279n/axMXFsWfPHvz8/Ni6dSvjxo1j7NixVKhQgTfeeIM333yTqlWr5vSvqeWx24XX19eXsWPHUqJECQ4fPkzDhg1ZvHgxSinatm3LtGnTWLlyJQkJCXh6elKnTh2WLFly5/3x8fF07dqVmJgYkpOT+eqrr+jatWuu5c6J4u4NfKiUWkbaslxxerxdyyuxsbH8/vvvLFq0iH379mFpaUnHjh2ZOnUqHTp0wNXVNdeO7ezszLPPPsuzz6YtyxodHc26detYtmwZEydO5KuvvqJFixa8+eab9OzZk6JFi+ZalsJg0KBBHDx4MEf36enpycyZM7Pc/sCBAxw5cgR3d3datmyJv78/rVq1urN98uTJzJ49+7457ezs+PPPPylatCiXL1+mWbNmdOnSJdf+ysvKVMilwC6ghlIqQinVRynVXynVP73JOuAkEE7a+pnv50pSTcsgMjKS4cOHU758eT744AOSkpKYPn0658+f5++//6ZXr165Wtjvp2TJkrz55pusX7+ec+fO8fXXXxMbG0u/fv2oUKECn376KVFRUXmaSctZTZo0oVy5clhYWODp6cnp06ez/F4R4ZNPPqF+/fq0b9+e8+fPc+nSpVzLmpXZMv/LZLsAH+RYIk17iLCwMKZMmcKvv/5KSkoKPXv2ZOjQoXh5mdflFe7u7owYMYLhw4ezZ88epk2bxqRJk5g+fTp9+vRh2LBh6FteP5pH6WHnFltb2zs/W1pakpKSkuX3LlmyhOjoaAIDA7G2tqZSpUq5euWtvreMli9ERkby9ttvU6NGDX777Tf69u1LWFgYS5YsMbvCnpFSimbNmrFy5UpCQ0Pp3bs3P/74I1WrVmXAgAFER0cbHVHLYdbW1iQnJ//n9bi4OEqVKoW1tTU+Pj6cOZOlO/dmmy7umlm7efMm48ePp3r16vz+++8MGTKEM2fO8P3331OlShWj4z2SGjVqMH/+fE6dOsWAAQP46aefqFq1KlOnTuXWrVtGx9NyyHvvvUf9+vXp3bv3Xa/37t2bgIAAGjVqxJIlS6hZs2buBhERQx4NGzYUTXsQk8kkS5YskfLlywsgL7/8soSHhxsdK0eFhobKCy+8IIBUqVJFVq1aJSaTyehYZiUkJMToCIa63+8PBEgWaqzuuWtm58yZM3To0IHevXtTsmRJ/Pz8WLlyJU888YTR0XJUzZo1Wbt2LZs2baJIkSK8/PLLdOvWjQsX9GQz7fHp4q6ZDZPJxNy5c6lbty47d+5kzpw57Nu3j9atWxsdLVc988wzHDhwgGnTprFp0yZq167NL7/88sCrZDUtK3Rx18zCqVOnaN++PQMGDKBZs2YcOnSIAQMGYGFROP4XtbKyYujQoQQFBVG3bl3eeustOnXqREREhNHRtHyqcPzL0cza0qVLqV+/PgEBAfz4449s2rSp0E4TrF69On5+fnz77bf4+PhQt25d1qxZY3QsLR/SxV0zTGJiIgMGDKBXr154eHhw+PBh3n333UJ/XxYLCws+/vhjgoODqVq1Ki+++CJDhw697/Q6TXsQXdw1Q4SHh9OiRQvmzp3LiBEj8PHx0XcKvUfVqlXx9/fnww8/ZPr06bRu3ZqzZ88aHUvLJ3Rx1/Lc6tWradiwIadPn+bvv//m66+/1vdIfwBbW1u+++47VqxYwZEjR2jQoAHr1683OpaWD+jiruUZEWHChAm8/PLL1KpViwMHDtCpUyejY+ULPXr0IDAwkPLly9OpUydmzpypZ9NoD6WLu5Ynbt26xZtvvslnn33Ga6+9hp+fHxUrVjQ6Vr5SrVo1/P396datG4MHD+b999/X4/B5ZN68eZQpUwZPT088PDzo0aMHp06deuh7Vq1aRdOmTfHw8KBRo0Zs3Lgxj9Kmy8qVTrnx0FeoFh7R0dHSqlUrAWT8+PH6KszHlJqaKiNHjhRAnnnmGYmJiTE6Uq4xlytU33//ffnxxx/vPP/hhx/E09Pzgf8vL1myRJo1ayYXLlwQEZHjx4+Lm5ubnDt37pGOq69Q1czW0aNHadq0KQEBASxfvpzPPvus0M+GeVwWFhZMnjyZhQsX4uPjQ4sWLTLtRWqP59ChQ9St+/8L0fXv35+LFy9y7ty5/7S9ceMGo0aNYsWKFZQuXRpI+6urbdu2bNmyJc8y58RiHZp2XwEBAXTo0AFLS0t8fX1p2rSp0ZEKlLfffpvKlSvz0ksv0bJlS/7991/q1KljdKxcM2gQ5PBaHXh6QlbuJHz48OH/fLb29vbExMT8Z5bXsmXL8PLyonz58ne9bmtry82bNx87c1bpnruWK/z8/Hj66adxcnJi586durDnkrZt27J9+3YAWrduzd69ew1OVPCcO3cOJyenu1bSSk5OvnMPoD59+tC9e/c72w4fPoyHh8d/9hMUFETNmjUJDQ2lf//+dO/enR9++CHXcuueu5bj/vnnH7p3707lypX5999/KVu2rNGRCrQ6deqwY8cOnnnmGdq1a4e3tzdPPfWU0bFynFFrdQQHB/+n1/7zzz/z9NNP4+HhwYIFC+4q7kWLFiUpKemu9rt27eLGjRu0adMGCwsL5s6di8lk4t1338213LrnruWopUuX0q1bN+rWrcu2bdt0Yc8jVapUYceOHVSsWJGOHTvi7e1tdKQC497x9k2bNjFp0iSmTZt23/YvvPACK1asuLMQy/Hjx+nbty8LFy68c68kb29vWrVqRbt27XIveFbOuubGQ8+WKXgWLlwoSilp06aNxMXFGR2nULp8+bI0adJELC0tZfny5UbHeWzmMFumV69e4u7uLl5eXtKgQQPp1auXHD169K42L7/88l3PFy5cKPXr15caNWpI1apVZefOnffd9/PPP//QYz/ObBld3LUc8euvv4pSSp599lm5efOm0XEKtWvXrkmrVq3E0tJSVq9ebXScx2IOxf1hLl++LP369ZMqVarIxIkT/7P94sWLUrt2bQkKCrrzmo+Pj3z00Ufy3nvvyezZsx+6/8cp7koMusqtUaNGEhAQYMixtZy1bNkyevfuTdu2bVm7di329vZGRyr0rl+/znPPPUdAQACrVq2ic+fORkfKltDQUGrVqmV0DMPc7/dXSgWKSKPM3qvH3LXHsmrVKl577TVatWqFt7e3LuxmwsnJifXr1+Ph4UH37t3ZsGGD0ZG0PKaLu5Zt3t7e9OzZk6ZNm7J27VocHByMjqRl4OzsfGdlp27durF582ajI2l5SBd3LVs2b95Mjx498PLyYt26dTg5ORkdSbsPFxcX/v33X6pXr06XLl3YuXOn0ZG0PKKLu/bI9u/fz4svvkiNGjXYsGEDzs7ORkfSHqJEiRJs3ryZsmXL0rlzZ0JDQ42OpOUBXdy1R3LixAk6duyIq6srGzZswMXFxehIWhaUKlWKjRs3Ym1tzXPPPafXZi0EdHHXsuzSpUs899xzpKSksHHjRtzd3Y2OpD2CKlWqsH79emJjY+nQoQMxMTFGR9JykS7uWpZcv36d559/nsjISP755x9q1qxpdCQtGxo0aMCaNWs4fvw4Xbp0ISEhwehIWi7RxV3LVHJyMi+//DJBQUH88ccfNGvWzOhI2mN4+umn+e233/D396dXr16kpqYaHalQ8vX1zdUT3Lq4aw8lInz00Uf8+++//Pjjj7zwwgtGR9JywKuvvsr06dNZs2YNo0ePNjpOoaSLu2ao7777jnnz5jFixAjeeecdo+NoOWjgwIH079+fqVOn8vPPPxsdx2wtXryYJk2a4OnpSb9+/Thz5gzVqlXj8uXLmEwmnnzySTZt2gRAt27daNiwIXXq1OHHH3+8s48NGzbg5eWFh4cH7dq14/Tp08ydO5cZM2bg6el557bNOSlLtx9QSnUAvgUsgfkiMvme7RWAX4Bi6W1Gici6h+1T337A/K1fv55OnTrRuXNnVq9efeeOdlrBkZycTMeOHdm2bRubN2+mdevWRke6y12X3xuwWkdoaCgjRoxg9erVWFtb8/7779OsWTOSkpLYsGEDTZs2JTw8nHnz5gFw9epVXF1dSUhIoHHjxvj5+WEymfDy8mLbtm1Urlz5TpuxY8fi6OjIsGHDHnr87N5+INP7uSulLIHvgWeACGCfUspbREIyNPsMWCEiPyilagPrgEqZ7VszX0eOHOHVV1+lfv36LF68WBf2Asra2po//viD5s2b89JLL7Fnzx6eeOIJo2OZjS1bthAYGEjjxo0BSEhIoFSpUowdO5Y//viDuXPncjDDF86sWbP4888/gbRFPsLCwoiOjqZ169ZUrlwZAFdX1zzJnpXFOpoA4SJyEkAptQzoCmQs7gLcXqbEGYjMyZBa3oqOjqZz5844ODjw999/4+joaHQkLRe5uLiwdu1amjZtSufOndm1a5d5XphmwGodIsKbb77JpEmT7nr95s2bd64ViI+Px8nJCV9fXzZv3syuXbsoUqQIbdu2JTExMe0OjQasG5yV7lhZIOMqsBHpr2U0FnhNKRVBWq/9o/vtSCn1nlIqQCkVcPtG9pp5uT0z5sKFC3h7e1OuXDmjI2l5oGrVqqxatYqwsDB69uypZ9Cka9euHStXriQqKgpIG3Y5c+YMI0eOpHfv3nz55Zd3VlOKi4vDxcWFIkWKcPToUXbv3g1A8+bN8fPzu7OI+dWrV4G0m7tdv34998Jndk9goAdp4+y3n78OfHdPmyHA0PSfm5PWq7d42H71/dzN08CBAwWQJUuWGB1FM8DcuXMFkDFjxhgdRUTM437uy5YtEw8PD6lXr554eXmJr6+vNG3aVFJSUkRE5MUXX5SFCxdKYmKidOjQQerVqyfdu3eXNm3aiI+Pj4iIrFu3Tjw9PaV+/frSvn17ERE5duyY1KtXTzw8PGTbtm33PXauLtaRXqw3Zng+Ghh9T5sjQPkMz08CpR62X13czc+SJUsEkEGDBhkdRTOIyWSSd955RwD566+/jI5jFsXdSI9T3LMyLLMPqKaUqqyUsgF6Avcu0HgWaAeglKoF2AF63CUfCQ4Opm/fvrRu3ZopU6YYHUcziFKK77//noYNG/L6668TFhZmdCQtmzIt7iKSAnwIbARCSZsVc0Qp9aVSqkt6s6HAu0qpIGAp8Fb6N4yWD8TExPDSSy/h4uLC8uXLsba2NjqSZiA7OztWrVqFtbU1L730EvHx8UZH0rIhK7NlkLQ56+vuee3zDD+HAC1zNpqWF0wmE6+99hpnz57Fz8+P0qVLGx1JMwMVK1Zk6dKldOjQgb59+7J06VJDZnwAhs02Mdrj9o/15OVCbsKECaxbt46ZM2fSvHlzo+NoZuSZZ55hwoQJLF++nG+//daQDHZ2dly5cuWxC11+IyJcuXIFOzu7bO9DL5BdiPn6+tKuXTt69erFr7/+Wih7R9rDiQgvvvgi69atw9/f/87FPHklOTmZiIgIEhMT8/S45sDOzo5y5cr9Z5g0q1eo6uJeSEVHR+Ph4YGTkxOBgYH6QiXtga5evUqDBg2wtLTkwIED5nmBUyGS1eKuh2UKIZPJxBtvvMHVq1dZsWKFLuzaQ7m6urJs2TLOnj1L3759C90QSX6li3shNG3aNDZs2MDMmTPx8PAwOo6WDzRv3pyJEyeycuXKOzfJ0sybHpYpZHbu3Enr1q156aWXWL58uR5n17LMZDLxwgsv4OPjw+7du/H09DQ6UqGkx9y1/9Bjp9rjio6OxtPTE0dHRwICAnBycjI6UqGjx9y1u4gI/fr1IzIykuXLl+vCrmVLyZIlWbJkCWFhYQwaNMjoONpD6OJeSPz222+sXLmS8ePH5/l0Nq1gadu2LaNGjWLhwoWsXr3a6DjaA+hhmULg1KlTeHh40KBBA7Zu3YqlpaXRkbR8LikpiRYtWnDq1CkOHTqEu7u70ZEKDT0sowGQkpLC66+/jlKKX3/9VRd2LUfY2NiwePFiEhISePvttzGZTEZH0u6hi3sB9/XXX+Pv78+cOXOoWLGi0XG0AqRmzZp88803bNq0idmzZxsdR7uHHpYpwPbt20eLFi3o0aMHv//+u9FxtAJIROjcuTObN28mMDCQOnXqGB2pwNNTIQu5Gzdu0KBBAxITEwkKCsLFxcXoSFoBdenSJerVq4e7uzt79uzB1tbW6EgFmh5zL+Q++eQTwsLC+OWXX3Rh13KVm5sbCxYsICgoiPHjxxsdR0uni3sB5Ofnx6xZs/j444956qmnjI6jFQKdO3fmrbfeYvLkyei/yM2DHpYpYOLj46lfvz4WFhYEBQXh4OBgdCStkIiNjaVu3bo4Ozuzf/9+PTyTS/SwTCE1cuRITp8+zc8//6wLu5anihUrxvz58wkJCWHs2LFGxyn0dHEvQLZu3cqcOXMYNGgQTz75pNFxtEKoQ4cO9OnThylTprBnzx6j4xRqelimgLh+/Tr16tXDxsaGgwcPUqRIEaMjaYVUXFwc9erVw8HBgf3792Nvb290pAJFD8sUMsOHD+fs2bMsWrRIF3bNUM7OzixYsICjR48yZswYo+MUWrq4FwA+Pj7MmzePIUOG0KJFC6PjaBrPPPMM7733HjNmzGDv3r1GxymU9LBMPpeQkED9+vUREYKDg3WvXTMbcXFx1KlTB1dXVwICArCxsTE6UoGgh2UKiXHjxhEeHs6PP/6oC7tmVpydnZkzZw6HDh1i6tSpRscpdHRxz8f279/PtGnT6NOnD08//bTRcTTtP7p06cIrr7zCl19+ydGjR42OU6joYZl8KiUlhSZNmnDhwgVCQkL0LQY0s3Xp0iVq1apFnTp18PPzw8JC9ykfhx6WKeCmT5/OgQMHmD17ti7smllzc3Nj+vTp7Nixg3nz5hkdp9DQPfd8KCwsjPr169OhQwdWr16NUsroSJr2UCLCs88+y549ewgJCaFcuXJGR8q3dM+9gBIR+vfvj42NDd9//70u7Fq+oJRi3rx5pKSk8MEHHxgdp1DQxT2f+f3339m6dSuTJ0/W61Zq+UqVKlUYN24c3t7e/PXXX0bHKfD0sEw+EhsbS82aNalQoQK7du3S66Fq+U5ycjJeXl7ExcUREhKCo6Oj0ZHynRwdllFKdVBKHVNKhSulRj2gzStKqRCl1BGllF7TLRd89tlnREdHM3fuXF3YtXzJ2tqaefPmce7cOcaNG2d0nAIt0+KulLIEvgc6ArWB/ymlat/TphowGmgpInWAQbmQtVALCAhgzpw5fPDBB3h5eRkdR9OyrUWLFrz77rvMmDGD4OBgo+MUWFnpuTcBwkXkpIgkAcuArve0eRf4XkRiAEQkKmdjFm6pqakMGDAANzc3vYyZViBMnjwZFxcX+vfvj8lkMjpOgZSV4l4WOJfheUT6axlVB6orpfyVUruVUh3utyOl1HtKqQClVEB0dHT2EhdC8+bNIyAggOnTp+Ps7Gx0HE17bK6urnzzzTfs2rWLBQsWGB2nQMr0hKpSqgfwnIj0TX/+OtBERD7K0GYtkAy8ApQDtgN1RST2QfvVJ1Sz5uLFi9SsWZNGjRrx77//6qmPWoEhIjz11FMEBwdz9OhRSpUqZXSkfCEnT6hGAOUzPC8HRN6nzV8ikiwip4BjQLWshtUebPjw4SQkJDBnzhxd2LUCRSnFDz/8QHx8PCNHjjQ6ToGTleK+D6imlKqslLIBegLe97RZAzwFoJQqQdowzcmcDFoY+fj4sHjxYkaOHEn16tWNjqNpOa5WrVoMGTKERYsWsXv3bqPjFChZmueulHoemAlYAgtFZIJS6ksgQES8VVqX8hugA5AKTBCRZQ/bpx6WebikpCQ8PDxISkri8OHDeqkyrcCKj4+nRo0alClThj179uhpvpnI0XnuIrJORKqLyBMiMiH9tc9FxDv9ZxGRISJSW0TqZVbYtcx98803HD16lNmzZ+vCrhVojo6OTJ06lcDAQBYuXGh0nAJDX6Fqhs6ePUvNmjV5/vnnWblypdFxNC3XiQht2rQhJCSE48eQOQJQAAAgAElEQVSP4+rqanQks6VvHJaPjRgxAki7ra+mFQZKKb777jtiYmL4/PPPjY5TIOjibma2bdvG8uXLGTlyJBUqVDA6jqblGQ8PDwYMGMAPP/xAUFCQ0XHyPT0sY0ZSU1Np2LAhMTExhIaG6jVRtULn6tWr1KhRg1q1auHn56en/96HHpbJh+bPn09QUBDTpk3ThV0rlFxdXZk4cSLbt29n6dKlRsfJ13TP3UzExMRQrVo16tati4+Pj+6xaIVWamoqzZo1IzIykqNHj+Lk5GR0JLOie+75zBdffEFMTAyzZs3ShV0r1CwtLZk9ezaRkZF89dVXRsfJt3RxNwOHDx9mzpw59OvXj/r16xsdR9MM17RpU9566y1mzJjBsWPHjI6TL+nibjARYdCgQRQtWlTfzlfTMpg8eTL29vYMHjzY6Cj5ki7uBluzZg1btmzhyy+/pHjx4kbH0TSz4ebmxueff8769evZuHGj0XHyHX1C1UCJiYnUrl0bBwcHDhw4gJWVldGRNM2s3Lp1izp16mBra0tQUJD+N4I+oZovfPPNN5w6dYqZM2fq/2k17T5sbW2ZMmUKISEhzJ8/3+g4+YruuRvk/PnzVK9eneeee47Vq1cbHUfTzNbtRT1CQkIICwsr9KuR6Z67mRszZgwpKSlMmzbN6CiaZtaUUkyfPp3Lly8zYcIEo+PkG7q4G+DgwYMsWrSIjz76iCpVqhgdR9PMnpeXF2+88QbffvstJ0/qdYCyQhf3PCYiDBs2DBcXFz799FOj42havjFx4kSsrKz0knxZpIt7Hlu3bh1btmzhiy++wMXFxeg4mpZvuLu7M3LkSFauXMn27duNjmP29AnVPJSSkkL9+vVJSUnh8OHD2NjYGB1J0/KVGzdu3LUkn4VF4euf6hOqZmj+/PmEhoby9ddf68Kuadng4ODApEmTCAgIYMmSJUbHMWu6555Hrl27RtWqValZs6a+T7WmPQaTyUTTpk25cOECx44dw8HBwehIeUr33M3M5MmTiY6OZvr06bqwa9pjsLCwYMaMGZw/f15PJX4IXdzzwNmzZ5kxYwa9e/emUaNMv3A1TctEq1at6N69O1OmTOHChQtGxzFLurjngU8//RQR0RdgaFoOmjRpEklJSYwbN87oKGZJF/dcFhAQwOLFixk8eDAVK1Y0Oo6mFRhVq1alf//+zJ8/X9/z/T70CdVcJCK0bduW0NBQwsPDKVq0qNGRNK1AiYqKomrVqrRv377Q3KNJn1A1A3///Tfbtm1j3LhxurBrWi4oVaoUI0aM4M8//8Tf39/oOGZF99xzSUpKCh4eHncuWLK2tjY6kqYVSDdu3KBatWpUrlyZHTt2FPjZaLrnbrBff/2VkJAQJk6cqAu7puUiBwcHxo0bx86dO/nrr7+MjmM2dM89FyQkJFCtWjXKli3L7t27C3xPQtOMlpKSQr169RARDh8+XKAXv9E9dwPNnj2b8+fPM2XKFF3YNS0PWFlZMXnyZI4dO8aCBQuMjmMWslTclVIdlFLHlFLhSqlRD2nXXSklSqlCe6VOTEwMEydOpGPHjrRp08boOJpWaHTp0oWWLVsyduxYbty4YXQcw2Va3JVSlsD3QEegNvA/pVTt+7RzAj4G9uR0yPxk8uTJxMXFMWnSJKOjaFqhopRi6tSpXLx4kenTpxsdx3BZ6bk3AcJF5KSIJAHLgK73aTcemAIk5mC+fCUiIoJZs2bx2muv4eHhYXQcTSt0mjdvzksvvcSUKVOIiooyOo6hslLcywLnMjyPSH/tDqVUA6C8iKx92I6UUu8ppQKUUgHR0dGPHNbcjR07FpPJxJdffml0FE0rtCZOnEhCQkKh/3eYleJ+vzOCd6bYKKUsgBnA0Mx2JCI/ikgjEWlUsmTJrKfMB0JCQvj55595//33qVSpktFxNK3QqlGjBu+99x7z5s0jPDzc6DiGyUpxjwDKZ3heDojM8NwJqAv4KqVOA80A78J2UvWTTz7BwcFBr4uqaWZgzJgxWFtb88UXXxgdxTBZKe77gGpKqcpKKRugJ+B9e6OIxIlICRGpJCKVgN1AFxEpmJPY7+P2xRMjRoygRIkSRsfRtEKvTJkyDBw4kKVLlxIcHGx0HENkWtxFJAX4ENgIhAIrROSIUupLpVSX3A5o7kSEkSNH4ubmxuDBg42Oo2lauhEjRuDs7Fxo/5rO0mVcIrIOWHfPa58/oG3bx4+Vf6xdu5YdO3bwww8/FLrlvszZzah4ooMiiQ2J5GZ4JEmnI+HqVSziYrCMj8XmRgx2iXFYpSZinXoLa9MtbEyJWJKCCQsEC1KVJaIsSFE2JFo5csvGkWQbR5LtHElxdMFUvCTKrRTW7iWxr1iKYnXLUbpJBexc7I3+9TXAxcWFESNG8Mknn+Dv70/Lli2NjpSn9O0HHoPJZMLT05OEhARCQkL0PWTyWPyF65z55zCxu4+ScuwENmfDcb5yAveb4RST2P+0T8GSWOXCdSsXbloXI9HOmRRre0xWtqRa2yLWtpisrEEEZUoFkwllMmGRnIjVrRvY3IrHJjke++TrOKVcxVWuYInpP8eJtihFlH1FrrlU4lalGtjUr4lri5qUb18DBzfHvPhotHQ3btzgiSeeoEaNGvj6+haIK8azevuBgnsDhjywfPlyDh06xJIlS3Rhz00iXAyI4MzKfSTuCMD+xCHKXDlM+ZTT1ElvkoIlEVaVuOz8BEEV/4epQiWsKrhTpKo7zrXcKV6vDMXKF6WEhSKnzoqkJqVy+WQMscejuH4iihtHz5EcfgaL82dxiD5D6UsHKR+xGqsdqTAn7T0RlhU4X7IBibUaUKSVF+W7NMDNqyzKIv8XHXPk4ODAmDFj+PDDD9m4cSMdOnQwOlKe0T33bEpJSaF27drY2toSFBSEhYW+TU9OSbp+i2O/7eXqX9uxD95Dxai9uJkuApCMFadtaxBVqh5J1etSpGk9SrWpRdmWlbBxML8v2FvXkzjrc4Lo7aHcOngUm+OHKX3xAJWTjmGRPqP4kkVpTpVpwa2GLSjRpSXVe3ph7WBjcPKCIykpiZo1a1KsWDECAgLy/b/VrPbcdXHPpgULFtC3b1/+/PNPunXrZnScfC35RhKhC3cRs9oH54N+1IjdjX36hc4nrWtwvlwTUr2aULxjE6p198DO2dbgxI/v+oV4Tv4ZRMzW/VgF7qVChD8VUk4BkIgtocVacK1xO0r+rz01ejXE0lb/kf04fvvtN9544w2WL1/OK6+8YnScx6KLey66desW1atXx83NjT179hSIcby8FrHjNKd+2ICNzwZqX9iCE/GYUByz9+RSjTbYPdeGau88SfHqxY2OmmcuHrjAqSU7ubXVn9KhPtRMPAhAHM4cLfMUSe06Um3QC5RuWDaTPWn3Sk1NxcPDg+TkZI4cOZKvbwmc1eKOiBjyaNiwoeRXs2bNEkA2bdpkdJR8w5RqkpBf9opvi9ESbltLBERAzlpWFN9a/WXniD8l5lSM0THNStSRKNnx8XLxrfGunLWsdOczCyniJT6tP5eQX/eJKdVkdMx8Y82aNQLI/PnzjY7yWIAAyUKN1T33R3T77HvNmjXx8fHRvfaHMKWYOPzDdmLnr6TqkTW4p0aQgiXBxdoQ16YzFd7rSJUO1fXJxCwQkxDuHULE3LW47vybutd3YYmJCMuKnGjQnZLv96DWm030Z/kQIkLz5s05f/48YWFh2NnZGR0pW3TPPZdMnjxZANmxY4fRUcxW+JpD4ttspERYlhcBuYmd7CrdVbb1WSSXj102Ol6BEB0aLdve+Vn2lHxebmEtAhJhWUF8Gw2VsNXBRsczW1u2bBFApk+fbnSUbEP33HNebGwsVapUoVmzZqxbty7zNxQiFwMiODZ2KWW2LKZ6YjApWHKgxLPc6t4bzy+64lhaz+/OLbGnYjg0wRtb7z9oEL0Ra1IILeJFdMc3qTuxF67V9S0xMnrmmWc4ePAgJ06coGjRokbHeWS6554LxowZI4AEBgYaHcUsJN1Ikl3D/pBA1/aSihIBCXZoKj4vz5JLhy4ZHa9QijoSJT4vfish9g1EQG5hLbvdX5Q9n6yRpBtJRsczC3v37hVAvvzyS6OjZAtZ7Lnr4p5FUVFR4ujoKC+//LLRUQwX4X9afFt9KhctSouAnLOsIFtbfyEnNx43OpqWwdEVQeLjNUQuKTcRkGhVUnxbjJbI3WeMjma4Ll26iLOzs1y9etXoKI9MF/ccNnToULGwsJAjR44YHcUQKbdSZO8Yb9lb6nlJRUkKFrK7VGfZ88U/knIrxeh42kMk3UyW3WPWyi63rpKCRdp/O/cXZf+0LYV2ts3BgwcFkM8++8zoKI9MF/ccdP78ebGzs5M33njD6Ch57lpEnPh2+UbOWVYQAblgUUa2thoj5/x17y8/Orv9tGxtOkouq+IiIOE2tcTv1e/l2vlrRkfLcz169BBHR0eJjo42Osoj0cU9Bw0YMECsrKzkxIkTRkfJM5F7zopv42ESS1ERkP3ObcR/2Co9bltA3LyaINv6LJLDRRqJgMThJL6NhsqFgAijo+WZI0eOiFJKRowYYXSUR6KLew45ceKEWFlZSf/+/Y2OkieOLt0vOyr1liSsJBlL2VGhpxxetM/oWFouMZlEgufvkR0Ve0kylnILa9lW/R05uS7U6Gh5onfv3mJvby8XLlwwOkqW6eKeQ9544w2xs7OTiIgC3KMxmeTA9K2y3+VpEZBrOIpPg8Fydvtpo5NpeeiM70nxqfeh3MBeUlGyu0w3OTx/l9GxctXx48fF0tJSBg4caHSULNPFPQeEhISIhYWFDBkyxOgoucKUapLAr/+Vg0WflNvj6T7PT5HY0/o2AIVZ1JEo2dpqjFxVLiIgB5zbSODX/xbYk69vv/222Nrayrlz54yOkiW6uOeA7t27i6Ojo0RFRRkdJUeZUk0SMHGjBDm1EAGJtCgrvt2/k5tXE4yOppmRa5HXxafLdDlvWS69yLeWAzN9jY6V406ePClWVlYyYMAAo6NkiS7uj+nAgQMCyJgxY4yOkmNMqSbZ9+U6OeTYTNIuVy8vvq/OkYTYRKOjaWYsITZRfHvMlgsWZURAAl2eluAfCtbtN/r16yfW1tZy+rT5D0Xq4v6YunbtKs7OzhITUzCGKIK+3y5BTi3l9p0Y/XrPk8Q4XdS1rLt55ab4dJ0hUaqUCEhA8Wfl8ILdRsfKEWfPnhUbGxvp06eP0VEypYv7YwgICBBAxo0bZ3SUx3Z89SHZ49ZJbo+p+/aaK7eu3zI6lpaPxV+KF58Xpkq0KiECssets4R75/+L+z766COxtLSUsLAwo6M8lC7uj6FTp07i4uIisbGxRkfJtgj/07LjiTckFSUxOMvWZydJfNQNo2NpBci1yOuytf0EiaWopGAh22q+K5cORhodK9siIyPF3t5eXn/9daOjPJQu7tl0+6ZCX331ldFRsuXK8cvi6zVYErGRBGzFp9EwuRJ2xehYWgEWHRotPp4D5RbWEk8R8Wn9uVy/cN3oWNkybNgwUUpJSEiI0VEeSBf3bOrYsaO4urrKtWv563LspBtJ4vviTIlRxdJ6UdXekYhdZ42OpRUipzaHi3+5HiIglyzcxK/XD5KckGx0rEcSFRUlDg4O8sorrxgd5YGyWtzz9zLgOWz37t2sX7+e4cOH4+TkZHScLAsYv56zLvVp8+cgTrg05uTqIJ48voCyzcobHU0rRCq1e4IW51ZweP5uLjhWo/XvAzjjXI/ASZuMjpZlJUuWZODAgaxYsYLg4GCj4zyerHwD5MbDHHvuzz77rJQoUUKuX88ff1KeWBsie0t2FAE5aV1Ndo/5u8BeaKLlL6ZUk+watUZOWz0hArK7dFc545M/7s105coVKVq0qHTr1s3oKPeF7rk/Gn9/fzZt2sSIESNwdDTvVYNiT8Xg12AgFTrVo1r0Tnw7f0PZq4dp+mUnvYamZhaUhaLZpK64XT6Cz3OTqHNxM6Weqo1vq8+4EXXD6HgP5erqypAhQ1izZg2BgYFGx8m+rHwD5MbD3Hru7dq1k1KlSkl8fLzRUR4oNTlVtr+zUKJVCUnBQnxr9ZOoIwXr6lmtYDq/N0K2V3pNBOS8ZTnx/3ipWf+VGRsbKy4uLvL8888bHeU/0D33rNu+fTtbtmxh5MiRODg4GB3nvsJWBXO4eGtaLXyH8441CFu2nzYhcylZu6TR0TQtU+6Ny9Lq1G8Efb+DOJtStJj1P4Jd2hD+1xGjo92Xs7Mzw4cPZ926dezdu9foONmTlW+A3HiYU8/9qaeeEjc3N7lxw/zmgV87f018vQZLMpYSrUqI39s/S2pyqtGxNC3bUm6liN9rP8plVVySsBKfZqPkRrQZ/tu7dk1cXV3NrveO7rlnja+vLz4+PowePZoiRYoYHecOMQm7Bq8gvnxNntw/E/9afbE4fozWC9/CwqrQ/2fT8jFLG0ta//YuhB5ld9XXaLt7MlfK1CFg3D9GR7uLk5MTw4YNy7+996x8AwAdgGNAODDqPtuHACFAMLAFqJjZPs2h524ymaR169ZSpkwZuXnzptFx7jjre0ICij8rAhJi30AOzS8Y9+/QtPs5MNNXwmxqiYDsKvuyRO4zn7UTbvfeX3jhBaOj3EFOXcQEWAIngCqADRAE1L6nzVNAkfSfBwDLM9uvORT3zZs3CyDfffed0VFERCQ5IVl8O0+TG9hLHE7i8/IsSU7Ui09rBV/itVvi0/4ruYmdXMNRfF+caTYLr0+YMEEA2bt3r9FRRCRni3tzYGOG56OB0Q9p3wDwz2y/Rhd3k8kkLVu2lLJly0pCgvH3MT+67IAcKdJQBGS3W2c5vyd/LBygaTnp9JZw2Vf8ORGQww5NJOzPQ0ZHkri4OLPqvWe1uGdl8LYscC7D84j01x6kD7D+fhuUUu8ppQKUUgHR0dFZOHTu+ffff/H39+fTTz/Fzs7OsBwJVxPwbT6aJ3o2okRiBDsHraBJ5F+4NylnWCZNM0rFp5+gYdR6/N9fgtvNk1R40Qvfp8aRFJ9kWKaiRYsydOhQ/vnnH/bt22dYjkeWWfUHegDzMzx/HfjuAW1fA3YDtpnt18ieu8lkkmbNmkn58uUlMdG4e5ofmOEjp6yrioD4VesjV09cNSyLppmb6JAo2V6xlwjIcdu6cmSRccMi5tR7Jwd77hFAxpuUlAMi722klGoPfAp0EZFbj/F9k+s2btzI7t27+eyzz7C1tc3z48dfjGdb3ffxHPwUSoT9U7fQ+vh8XKq45HkWTTNXJWqVpNXpJez51BvH5BhqvNUMv8bDuHn5Zp5nydh7DwgIyPPjZ0tm1R+wAk4Clfn/E6p17mnTgLSTrtWy8o0iBvbcTSaTNGnSRCpWrCi3buX9ohUHZvjIGavKkooSH68hZjm/V9PMTeyZWPGt1U8E5LT1E3Jwll+eZ4iLixMXFxfp1KlTnh87I3Kq5y4iKcCHwEYgFFghIkeUUl8qpbqkN5sKOAJ/KKUOKqW8c/QbKAdt2rSJvXv38umnn2JjY5Nnx70RdQM/z4/xHPwUJmXJodnbaBv4DUVKmM/cek0zV84VnGkTMpcD031AoN7HbfFrNISEqwl5luF2733t2rX5o/eelW+A3HgY0XM3mUzSvHlzqVChQp722oO+335nbN3H42OJv2S+96/RNHN3/WK8+NZ9XwQk3KZmno7Fm0PvHX2F6n9t3ryZXbt2MXr06DzptSdcTcCv0VDqftAaC0nlwAxf2h78FodS5nn/Gk3LDxzdHGhz6HsCJm7CPjWe6m81x/fJMXkyoyZf9d6z8g2QG4+87rnfntderly5PJkhc+inXXLSuroIiG+d9/PtsmOaZs5iT8fI9ifeFAEJtfeUYyuDc/2Yt3vvnTt3zvVj3Q+65343Hx8f/P39GTVqVK7OkEmKT8K3+ShqvdsSa1MigV9vps3h73Esbd73iNe0/Mi5YjFahS9i96g1FE+MpFL3hvh2nEzqrZRcO+bt3vvff/9t1vd7V2lfBHmvUaNGkpd/1rRp04bw8HBOnDiRaxctnfg7hORXe1Mz4SDbqvfFc8s3FC1XNFeOpWna3S4fvUxY+wE0P7+Sw47NcFz1C5WerZ4rx7p27RqVKlWiVatWeHvn7fwRpVSgiDTKrF2h6Ln7+fmxbds2Ro4cmSuFXUzCtldm496lISUSI9jzyV+0PvaTLuyalodK1CxBs7Mr2Pnh75S9cYxSz3myrfdcxJTzHdiiRYsyZMgQs+69F4qe+9NPP01oaCgnT57E3t4+R/cdFXSBs+3fodHlDewt+TyVNi+gVP3SOXoMTdMezcX9kUQ8+zaNrmxir1tnqvgsoEStnF3YJi4ujsqVK+d571333NNt374dHx8fRo4cmeOFfe/oP7FsUI86l33xe/V7Gl9cqwu7ppmB0l7ueF1cj1+3GXhc2oipbj0CJ2zI0WM4Ozubd+89K2ddc+ORV7Nl2rdvn+OrLF2/cF22Ve8jafdb95Lwv0NybN+apuWsY38EyXHbOmkz1zw/loSYnLsL7O21Vrt06ZJj+8wMerYM7Ny5k82bNzN8+PAcW2XpyILdXKngScvjC/FpPpononbxRKdaObJvTdNyXvXu9SkXuQ8/j49pc3AW50o35viqQzmy79u9d29vb/bv358j+8wxWfkGyI1HXvTcn3vuOSlZsqTExz/+FaEpt1LEp914ScZSzlpWlAOztuVAQk3T8tLecevkkoWbJGArvi/OzJH1iGNjY6VYsWJ51nunsPfc9+zZw8aNGxk2bBgODo93ReiFgPMEu7Wn7ZYx7KnwCkVPBuH50ZM5lFTTtLzS+POOqOBggtyepc2fgzhQuiOXDl54rH2abe89K98AufHI7Z57x44dpUSJEnL9+uNdGbpnjLdcVsXlOg6yrc8iMaWaciihpmlGMaWaxLfnD3IDe4lSJWXf+PWPtb+87L1TmHvu+/btY/369QwdOhRHx+xdGXorLhE/z4E0Gd+FaNvyRK0L5Mn5b6IsVA6n1TQtrykLRZul/Yn0DiTGpjSNxnTEt8mIbN+fxtnZmcGDB+Pt7c2BAwdyOG02ZeUbIDceudlz79Spk7i6usq1a9ey9f4T645KqL2nCIiP50BJiDVutSZN03LXzSs3xbf2ABGQQw5N5IzPiWztJyYmRpydneWll17K4YR3o7D23Pfv38/atWsZMmQITk5Oj/ReMQk7+i6i1PMNKZl4jj1j/qbtgZnYOef9ak2apuUNe1d72hyZw66hKyl/4xjOTzVg58Dlj7yfYsWKMXDgQFavXs2hQzkzG+exZOUbIDceudVz79q1qxQrVkxiY2Mf6X1x5+LEv+L/RED2O7eVyH0RuZJP0zTzdW77KQl2bC4C4lfz3UdeKe3KlSvi5OQkr7zySi4lLKQ994MHD/LXX38xePBgnJ2ds/y+kF/2EVO5AU3OrMCn3VfUj9pMmUZlczGppmnmqFyrStS85IdP89G0OjqfyHKNCfvzcJbf7+rqykcffcQff/xBSEhILibNXIEq7uPHj8fZ2ZmPP/44S+1NKSZ8X5hKtbdaYCkphMzx46nNn2JpY5nLSTVNM1fWRax5audEDkzaSNHkK5R7qTHbXvsxyzcgGzx4MEWKFOGrr77K5aQPV2CKe3BwMKtXr2bgwIEUK1Ys0/ZXjl0m0L0TbdeNIMC9K45hB6k/oGUeJNU0LT9oOOoZOBjEkeKtab2kH7srvkrcmdhM31eiRAk+/PBDli1bxtGjR/Mg6f0VmOI+fvx4ihYtyqBBgzJtG/yDP7fqNKB+9Bb8XvmeZuf+oFhllzxIqWlaflKqnhteF9fj2/FrGkX8ybWqDTi8YE+m7xsyZAj29vZMnDgxD1LeX4Eo7ocPH2blypV8/PHHuLg8uEibUkz4Pj+F2u+3IVnZcnLxLtosf1/PXdc07YEsrCxou24Ex37aDkD1vk/i9+LMhw7TlCpVigEDBrBkyRLCw8PzKurdsnLWNTceOTlb5tVXXxVHR0e5fPnyA9tcOX5Z9pR6QQRkZ9nuEnvm0WbTaJqmxZy8KrtLdxUB2VWmm8ScvPrAthcuXBA7Ozt5++23czQDhWW2TGhoKCtWrOCjjz6iePHi921z6MddJNRqgEfUv/j1mE2zsytwrpD12TSapmkAxSq70OT8n/h2nU7DC2u5Xt2LkF/23bdt6dKl6devH7/++isnT57M46QFYFjmq6++okiRIgwZMuQ/28Qk+HWeRs1+rUlVVpz4dSdtVnygh2E0Tcs2ZaFou2Ywx37ajhITVd9qiV/37+47TDNixAisrKyYNGlSnufM18X92LFjLFu2jA8++IASJUrctS3mxFX2unelzdrhBLh3xTl8P7Vfb2hQUk3TCpq6fZvhcOwAB9w60GbVx+yu0IO4s3F3tXF3d+fdd99l0aJFnDlzJk/z5eviPmHCBOzs7Bg6dOhdrx+ev5sbNRrQ4NIG/F6eRbNzf+BcMfPpkZqmaY/C5QlXmkT+he8LU2l8fg2xT3gRuuTu2/6OHDkSCwsLJk+enKfZ8m1xDwsLY8mSJQwYMIBSpUoBacMwvl2nU+PdJzEpS8IX+dNm5Ud6GEbTtFyjLBRt1w4jdO42rCWJKq81x6/nnDvDNOXKleOdd95hwYIFnDt3Ls9y5dviPmHCBGxtbRk+fDgAsSevsrdsN9p6DyWwTOe0YZg3GxucUtO0wqJevxbYhR4kuGR72iz/gF2VenIt4hoAo0aNQkSYMmVKnuXJl8X9xIkTLF68mP79++Pm5saRhXuIr+5Fg4vr8XtxJk0jVulhGE3T8pxrteI0jPwb3w6TaXJuFVeqNOLY8oNUrFiRt956i59++onIyMg8yZIvi/vEiROxtrZm2NBh+L04k2p90pa8C1u4gzarB+phGE3TDGNhZUHb9SMJ+d4Xu9SbVOzZjG295zF61GhSUlLyrPeu0ubE5xjEt6EAAAX9SURBVL1GjRpJQEDAI7/v1KlTVK9enQ97f0DPjadpevEv9pTuSo2dP+tbCGiaZlYuh0Zz+snXaXRlI/4V/8eiZkVY/NcSTp06RenSpbO1T6VUoIg0yqxdlnruSqkOSqljSqlwpdSo+2y3VUotT9++RylV6dEjZ83PP/9MY6nJ4MV/4XXxH3y7TqfJ+T91Ydc0zeyUqFUSr4vr8H1mAs3OLGf06m1UT6zAb7/9luvHzrS4K6Usge+BjkBt4H9Kqdr3NOsDxIhIVWAG8HVOB72tfUQlfFOPYYGJY/N30HbNYD0Mo2ma2bKwsqDtpk84NHMrRVLj2c1ZWkZWzP3jZqFNEyBcRE6KSBKwDOh6T5uuwC/pP68E2imlcqXiFmtSk/3unXE4doC6fZrmxiE0TdNynOfANlgGHyTY7RmKN6ue68ezykKbskDGyZkRwL1V9U4bEUlRSsUBxYHLGRsppd4D3gOoUKFCtgLX798C+rfI1ns1TdOMVLJOKUpe9M6TY2Wl536/Hvi9Z2Gz0gYR+VFEGolIo5IlS2Yln6ZpmpYNWSnuEUD5DM/LAfdO1LzTRillBTgDV3MioKZpmvboslLc9wHVlFKVlVI2QE/g3r8rvIE303/uDmwVo+ZYapqmaZmPuaePoX8IbAQsgYUickQp9SVpN433BhYAvymlwknrsffMzdCapmnaw2XlhCoisg5Yd89rn2f4ORHokbPRNE3TtOzKl7cf+L/27iZEjiqK4vj/iIgI8QMGQYIaAgYcZmPIIm78QBGZRdyIRAgiDEIiuHGbjehOUEEImCxEXCjqRgfBjZIQCY4SSExiQPEjiBBMFpqNKBGvi1eBYXS6n2PVe11V5wcN3dPF9D39qu/UVFXXMzOzydzczcwGyM3dzGyAql04TNJFYCPzTs2x5stRIzDGzDDO3GPMDOPMvdHMt0fE1C8KVWvuGyXpeM4V0YZkjJlhnLnHmBnGmbvrzN4tY2Y2QG7uZmYD1Mfmfqh2ARWMMTOMM/cYM8M4c3eauXf73M3MbLo+brmbmdkUbu5mZgM0s819luZtLSUj87OSzko6JekTSd3P1VXAtNyrlntUUkjq/SlzOZklPdaM91eS3ipdY9sy1u/bJB2WdKJZxxdr1NkmSa9LuiDpzDrPS9KrzXtyStL21l48ImbuRrr65HfAVuAa4Etgfs0yTwOvNfd3A+/UrrtA5vuB65r7+/qeOTd3s9wm4CiwAuyoXXeBsb4DOAHc1Dy+uXbdBTIfAvY19+eBc7XrbiH3PcB24Mw6zy8CH5EmPNoJfN7Wa8/qlvtMzdtayNTMEXE4In5rHq6QJk7pu5yxBngBeBH4vWRxHcnJ/BRwICJ+AYiIC4VrbFtO5gCub+7fwD8nBeqdiDjK5ImLHgHejGQFuFHSLW289qw293+bt3XzestExJ/AlXlb+yon82pLpL/4fTc1t6S7gFsj4sOShXUoZ6y3AdskHZO0IunhYtV1Iyfzc8AeST+RLjH+TJnSqvqvn/tsWddzr6C1eVt7JDuPpD3ADuDeTisqY2JuSVcBrwBPliqogJyxvpq0a+Y+0n9on0paiIhfO66tKzmZHwfeiIiXJN1NmgBoISL+6r68ajrrY7O65T7GeVtzMiPpQWA/sCsi/ihUW5em5d4ELABHJJ0j7Zdc7vlB1dz1+4OIuBwRPwBfk5p9X+VkXgLeBYiIz4BrSRfXGrKsz/1GzGpzH+O8rVMzN7snDpIae9/3wV4xMXdEXIqIuYjYEhFbSMcadkXE8TrltiJn/X6fdAAdSXOk3TTfF62yXTmZfwQeAJB0J6m5XyxaZXnLwBPNWTM7gUsRcb6V31z7aPKEo8yLwDekI+z7m589T/pgQxr494BvgS+ArbVrLpD5Y+Bn4GRzW65dc4nca5Y9Qs/PlskcawEvA2eB08Du2jUXyDwPHCOdSXMSeKh2zS1kfhs4D1wmbaUvAXuBvavG+UDznpxuc9325QfMzAZoVnfLmJnZ/+DmbmY2QG7uZmYD5OZuZjZAbu5mZgPk5m5mNkBu7mZmA/Q3w5JeUhxwVwgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import sympy as sp\n",
    "import pylab as plt\n",
    "import pylbm\n",
    "\n",
    "u, X, LA = sp.symbols('u, X, LA')\n",
    "\n",
    "def solution(x, t):\n",
    "    return np.sin(np.pi*x)*np.exp(-np.pi**2*mu*t)\n",
    "\n",
    "xmin, xmax = 0., 1.\n",
    "N = 128\n",
    "mu = 1.\n",
    "Tf = .1\n",
    "dx = (xmax-xmin)/N # spatial step\n",
    "la = 1./dx\n",
    "s1 = 2./(1+2*mu)\n",
    "s2 = 1.\n",
    "dico = {\n",
    "    'box':{'x': [xmin, xmax], 'label': 0},\n",
    "    'space_step': dx,\n",
    "    'scheme_velocity': la,\n",
    "    'schemes': [\n",
    "        {\n",
    "            'velocities': list(range(3)),\n",
    "            'conserved_moments': u,\n",
    "            'polynomials': [1, X/LA, X**2/(2*LA**2)],\n",
    "            'equilibrium': [u, 0., .5*u],\n",
    "            'relaxation_parameters': [0., s1, s2],\n",
    "        }\n",
    "    ],\n",
    "    'init': {u: (solution, (0.,))},\n",
    "    'boundary_conditions': {\n",
    "        0: {'method': {0: pylbm.bc.AntiBounceBack,}},\n",
    "    },\n",
    "    'parameters': {LA: la},\n",
    "    'generator': 'numpy'\n",
    "}\n",
    "\n",
    "sol = pylbm.Simulation(dico)\n",
    "x = sol.domain.x\n",
    "y = sol.m[u]\n",
    "\n",
    "plt.figure(1)\n",
    "plt.plot(x, y, 'k', label='initial')\n",
    "\n",
    "while sol.t < 0.1:\n",
    "    sol.one_time_step()\n",
    "    \n",
    "plt.plot(x, sol.m[u], 'b', label=r'$D_1Q_3$')\n",
    "plt.plot(x, solution(x, sol.t),'r', label='exact')\n",
    "plt.title('Heat equation t={0:5.3f}'.format(sol.t))\n",
    "plt.legend();"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
